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Abstract 

Mean field approximation of a large collection of FitzHugh-Nagumo 
excitable neurons with noise and all-to-all coupling with explicit time- 
delays, modelled hy N ^ 1 stochastic delay-differential equations is 
derived. The resulting approximation contains only two deterministic 
delay-differential equations but provides excellent predictions concern- 
ing the stability and bifurcations of the averaged global variables of 
the exact large system. 
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1 Introduction 



Small parts of brain cortex may contain thousands of morphologically and 
functionally similar interconnected neurons. Realistic models of an individual 
neuron, like Ho dgkin- Huxley, FitzHugh-Nagumo (FN) or Hindmarsh-Rose to 
mention only a few popular examples [1], are given by few-dimensional non- 
linear differential equations. Transport of information between neurons can 
be phenomenologically described by time-delayed inter-neuronal interaction 
(please see [2] and the references therein). It is also well known that neurons 
in vivo function under influences of many sources of noise [3j. Considering all 
mentioned factors it is clear that a basic, relatively detailed mathematical 
model of a small part of realistic cortex should involve an extremely large 
system of nonlinear stochastic delay-differential equations (SDDE). Analyzes 
of such complex models is impossible without more or less severe approxima- 
tions, which should be adopted to different purposes. It is our goal to study 
some aspects of an approximation by only two deterministic delay-differential 
equations (DDDE) of an example of a complex neuronal system described by 
many-component SDDE. We shall see that, although the approximate model 
is very simple, the predicted critical parameter values for the bifurcations and 
stability of the stationary states are in excellent quantitative agrement with 
those of the exact complex model within a relevant domain of parameters. 

Neuronal dynamics with all three factors (large number of units, delayed 
interaction and noisy environment ) included has been studied much less 
than the influence of each of the factors separately [1]. Important influence 
of noise alone on a single, small number or large clusters of neurons has 
been studied a lot in recent years [S]. It is also well known that time-delay 
can have important qualitative effects on the stability of stationary states 
(please see for example [2], [6]) and synchronization of neuronal dynamics 
[7] . Studies of combined effects of noise and time-delay have mostly, but not 
entirely ([H]) been restricted to artificial networks [8] [9] or small number of 
neurons (usually two) [IO],[TT]. An example of a study of a large collection 
of noisy realistic neurons with delayed coupling can be found in [T7] (see also 
the references therein). 

The mean field approach (MFA) is based on a set of approximations that 
replace many component system by a simpler system described by a small 
number of (averaged) collective or macroscopic properties. The mean field 
approximation has been applied on systems of excitable neurons with noise 
but with no time-delay for example in [12] , [13] , [S] , [S] • On the other hand 
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a type of MFA was devised in and [TB] and applied on large clusters of 
noisy neurons with time-delayed interaction in [T7]. However, the approxi- 
mations made in these papers resulted in a system of equations that is still 
to large to be analyzed analytically, so that the approximate system must 
be studied numerically. We shall derive an approximate system of only two 
DDDE for the dynamics of the mean fields. Such a simple system allows 
analytical treatment of bifurcations and the parameter domains of stability 
of the stationary states which turn out to be in a quite good agrement with 
the exact complex system. 

2 The model and its mean field approxima- 
tion 

We shall study a system of excitable neurons modelled by the following set 
of SDDE: 

c ^ 

edxi = f{xi,yi)dt + —"^{xjit - t) - Xi)dt 

i=i 

dVi = g{xi,yi) + V2DdW, (1) 

with 

f{x,y) = X - x^/3 - y + I 

g{x,y) = x + b, (2) 

where b, I, c, D and e ^ 1 are parameters. The formulas (2) represent one 
of the common ways of writing the famous FitzHugh-Nagumo model pQ of 
the excitable behavior. For certain parameter values, like b = 1.05, / = to 
be used throughout this paper, the ODE given by (2) have stable stationary 
solution {xo,yo) such that small departures from (xo,|/o) might lead to large 
and long lasting excursions away from {xo,yo) which nevertheless end up on 
the stable state {xo,yo). The type of excitable behavior epitomized by the 
FN model is called type II and is characterized by destabilization of the 
stationary state via the Hopf bifurcation. The variable x is called the fast 
variable (due to e ^ 1) and corresponds to the membrane electrical potential. 
The variable y is the slow recovery variable and has no direct interpretation. 

Each of i = 1, 2 . . . units in (1) is coupled with each other unit and with 
itself. There are two major types of inter-neuronal couplings: the chemical 
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and the electrical synapses. Time-delay r is important especially in the first 
type of synapses but plays also an important role in the electrical junctions 
and in the transmission of an impulse through the dendrite. In (1) we use 
the electrical coupling with the time-lag and the strength that is equal for 
all pairs of neurons. 

The terms \/2DdWi represent stochastic increments of independent Wiener 
processes, i.e. dWi satisfy 

E{dWi) = 0, E{dWidWj) = 6,jdt, (3) 

where E{) denotes the expectation over many realizations of the stochastic 
process. 

Mean field approximation 

In order to derive the approximate dynamical equations for the mean 
fields 

Xit) = -J2x,it)^<x.it)>, Y{t) = -Y.y,it)^<y,it)> (4) 

i i 

of the system (1) to be used in this paper we shall first suppose that: a) 
The dynamics is such that the distributions of Xi and yi are Gaussian and 
b) for large the average over of local random variables is given by the 
expectation with respect to the corresponding distribution, i.e. for example 
jfJ2fxi ~ E(xi), where E(xi) is the expectation with respect to the dis- 
tribution of Xi{t). In the limit N ^ oo the last assumption is expected to 
become an equality, implied by the strong low of large numbers |l8j. In the 
mean field approach it is commonly assumed that b) is approximately true 
even for finite but large N despite the nonzero interaction between the local 
random variables. The first assumption should be expected to be true when 
the noise intensity is small, i.e. D <^ 1 (see for example tl3j,[14j). With these 
assumptions the system (1) of 2N SDDE can be reduced to five DDDE for the 
macroscopic variables X{t),Y{t) and the second order cumulants. Further 
assumption concerning the time scales of first and second order cumulants 
enables us to derive the final approximate system of only two DDDE. 

Mean field assumption guaranties that global averages, like (1/A^) 
of local quantities are equal to the expectations with respect to distribution 
of the corresponding variable E{xi). Besides the mean values X{t),Y{t) we 
introduce deviations from the expectations: n^.it) = X{t) — Xi{t), ny.{t) = 
yi{t) — Y{t). Because of the assumed Gauss distribution of each variable 
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the first and the second order cumulants of these deviations are equal to the 
first and second order moments ( i.e. to the first and second order centered 
moments of the variables Xj, etc. . . ). Furthermore, due to the same Gaussian 
assumption higher order cumulants are equal to zero, and this enables us to 
terminate the cumulant expansion of the dynamical equations. Details of 
the derivation are given in the appendix. The result is a system of five 
deterministic delay-differential equations for the global variables and global 
centered moments: 

=< nl^it) >,Sy =< nl^{t) >,u =< n^Uy > . (5) 

The equations are 
dX{t) 



dt 
dYjt) 

dt 
e dsx(t) 

2 dt 

1 dSyit) 

2 dt 



X{t) - X{tf/?> - s^{t)X{t) - Y{t) + c{X{t - r) - X(t)), 
X{t) + 6, 

s.,{t){l - X{tf - s^{t) - c) - u{t) 
u{t) + D, 



^ = ^{l-X{tr-s.{t)-c)--^Sy{t) + s^{t). (6) 

The analogous set of ordinary differential equations was used to study 
the mean field approximation of the stochastic system of N neurons without 
delay in f5]. The equations (6) are delay-differential equations because the 
original system of stochastic eq. (1) contains time-delay. 

In order to further simplify the approximate system we shall suppose that 
relaxation time-scale of the second order moments is much faster then those 
of the first order moments. Thus we can replace in eq. (6) the stationary 
values of s^, Sy and u obtained by setting the right hand sides of the last 
three equations in (6) equal to zero. As the results we obtain the following 
two DDDE: 

= x{t)-X{tf/?,-^[l-c-X{tf + {{c-l + X{tff + ADY/^' 
- Y{t) + c{X{t-r)-X{t)), 

dY{t) 



dt + 



(7) 
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3 Stability and bifurcations of the stationary 
state 



Stationary states, their stability and local bifurcations of the approximate 
system of DDDE (7) are determined by the standard procedure. It is remark- 
able that such a crude approximation provides relevant information about the 
exact system. 

There is only one stationary state of (7) given by: 



X{t) = Xo = -b, Y{t) = Yo = — [l + 673 + c - (4D + (c + 6^ _ i)2)i/2 

(8) 

Local stability of (8) is determined from the roots of the characteristic equa- 
tion. Due to the time-delay the system (7) has an infinite-dimensional state 
space, and the characteristic equation is transcendental with an infinite num- 
ber of roots. The characteristic equation is 



1 + (m2 + 4D)i/2 



1 

Ye 
1 c 

A exp(— Ar) = 0. 



2b'^m 



(m2 + 4D)V2 



A 



(9) 



where m = c — 1 + 6^. 

Bifurcations of the stationary state occur for those values of the param- 
eters such that any of the infinite number of roots of (9) has the real part 
equal to zero [IS]. This is possible only if A = iu, where cu can be taken to 
be positive. Substitution of X = iu in (9) gives 



where 



e + c'/e' + 2/e ± {k' - c'/e' - 21 tf - A/t') 



c2U/2 



1/2 



(10) 



k 



2e 



26^m 



ADyn 



The critical values of the time-lag r are related to the other parameters c, D 
an b by 

cos-i(-A;e/c) + 2j7r)l j = 0, 1, 2 . . . (11) 



if 



-ul + 1/e 
cuj±/e 



> 



(12) 
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Figure 1: Bifurcation curves (r^j_,c) (fig. a,b,c) for fixed D = (a), D = 
O.OOf (b), D = 0.003 (c) and (<±,L') (fig. d,e,f) for fixed c = -0.05 (a) 
c = 0.05 (b) c = O.f . In all figures b = 1.05. Gray curves correspond to r^^. 
and black curves to Tc^_,_, for j = 0, 1, 2, 3, 4, 5. 

If (12) is not satisfied then 

- cos~\-ke/c) + (2j + 2)7rl /uj±, j = 0, 1, 2 . . . (13) 



It can be shown by direct substitution that 

so that on the bifurcation curves rl^j^ or Tc^_,_ one unstable direction is created 
or destroyed. Together with the stability properties for r = the bifurcation 

curves (11), (13) and (10) completely solve the problem of stability of the 
stationary state. It can be shown by rather lengthy calculations that the 
bifurcations for t^j^ are the Hopf supercritical or subcritical bifurcations of 
the DDDE (7). 

Bifurcation curves ± (c) for fixed D, 6 = 1.05 and <±(£>) for fixed c, 6 = 
1.05 are illustrated in figure 1 for different values of D (fig. la,b,c) and c 
(fig. ld,e,f). The value h = 1.05 renders the stationary state Xo,Yo stable 
and excitable when r = and D = 0. 
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c c D 



Figure 2: Enlarged parts of bifurcation diagrams presented in fig. la,c,e with 
parameter values, indicated by letters: a,b,c,d (fig. a),e,f (fig. c),g,li (fig. e), 
that are used for comparison with the exact system presented in fig. 3. 




Figure 3: Illustrates dynamics of the global variable X(t) for the exact sys- 
tem of N=95 units for parameter values corresponding to the stable or un- 
stable state of the approximate system (7). Parameter values correspond- 
ing to a,b,c,d,e,f,g,h, indicated in fig 2, are a) (c, r) = (—0.12,0.14); b) 
(-0.06,0.11), c) (-0.06,0.29), d)(-0.06, 0.59), e)(0.07, 0.09), f)(0.08, 0.27), 
g){D,T) = (0.002,0.02), h)(0.002,0.29). 
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The predictions given by the bifurcation values (11) and (13) of the system 
(7) are check versus the numerical solutions of the exact system (1). To 
check the approximate predictions of the bifurcations of stability for the 
noisy system, i.e. when D 7^ 0, a proper notion of stochastic bifurcations 
would be necessary [18]. Instead we use the sample paths of the SDDE 
(1) with large and for D ^ to illustrate that these paths remain in 
the vicinity of the stationary solution if the approximate system's stationary 
state is stable, or near a periodic solution when the state of the approximate 
system is unstable. Figure 2 presents enlarged parts of bifurcation diagrams 
in figure 1, where particular values of the parameters that correspond either 
to stable or to unstable stationary state of (7) are indicated. These parameter 
values are replaced in the original system (1) with large and particular 
sample paths of (1) with these parameter values are computed numerically. 
Time series of the global variable X{t) along such sample paths are shown 
in figure 3, for the system (1) with = 95. There is nothing special with 

= 95 and the same qualitative behavior of X{t),Y{t) is obtained for 
any moderately large A^. It is clear that when the bifurcation diagrams of 
the approximate system (7) predict that the stationary solution is stable ( 
like in the cases: b,d,f,h) the sample paths of the exact system display small 
stochastic fluctuations around the stationary state. On the other hand, when 
the stationary state of (7) is unstable, as shown in the bifurcation diagrams, 
the sample paths of the exact system displays coherent oscillations with large 
amplitude, indicating that the exact system has stochastically stable periodic 
solution. The quantitative agrement between the domains of stability in 
the parameters {D, c, r) space of the approximate system (7) and the exact 
system (1) is indeed quite remarkable. It should be expected that such an 
agrement should be observed for small values of D since this is one of the 
conditions which guaranties that the local random variables have Gaussian 
distribution which is one of the asumptions in the derivation of the mean field 
equations. It is interesting to observe the domains of the time-lag where the 
bifurcation diagrams in fig. 1 and fig. 2 predict that the non-zero time-lag 
induces stabilization of the stationary state. This secession of oscillations due 
to the specific non-zero interval of the time-lag values is correctly predicted 
for the global variables of the exact system. 

It should be stressed that the agrement in the predictions of the approxi- 
mate system and the large exact system goes only as far as the parameter do- 
mains of stability are considered. It should be expected that the predictions 
of the parameters stability domains based on (7) should well approximate 
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the parameters stability domains of the exact system for small values of D 
since this is one of the basic assumption in the derivation of the mean field 
equations. Also, interaction strength c should be relatively small in order for 
the mean field assumption to be valid for moderately large (but finite) N . 
However, this domain of small values of c includes interesting bifurcations 
predicted by (7) and occurring in (1). On the other hand, large values of 
r induce unstable stationary state and stable oscillatory behavior in both 
systems (7) and (1) so formally there is no restriction on the time-lag r. We 
should make clear that it should not be expected that the values of X and 
Y for the deterministic approximate system (7) should reproduce stochastic 
orbits X{t),Y{t) for the large exact system or their ensemble averages. The 
correspondence between the orbits of the two systems for the same values of 
the parameters in different domains is only qualitative in the sense that they 
share the same types of attractors. 

4 Summary 

We have studied validity of the mean field approximation for the treatment 
of stability and bifurcations of the stationary state of a large collection of 
FitzHugh-Nagumo excitable neurons with noise and all-to-all coupling with 
delays, modelled by ^ 1 stochastic delay-differential equations. Standard 
assumptions of the mean field approach are used to derive the system of only 
two deterministic delay-differential equations. The stability and bifurcations 
of the stationary state of the approximate system can be studied analytically. 
The bifurcation curves of the approximate system give relevant information 
about the global variables of the exact large system. For zero and sufficiently 
small noise there is remarkable quantitative agreement of the parameters 
bifurcation values. On the other hand, it should not be expected that the 
approximation gives applicable results when the noise is to large, primarily 
because the assumption about the Gaussian distribution of values of the 
dynamical variables is not valid for large noise. 

Using the approximate system it is predicted, and confirmed by direct 
numerical simulations on the large exact system, that the time-lag in a non- 
zero interval can stabilize the global variables onto the stationary values even 
when for zero time-lag the global variables perform large oscillations. This 
is reminiscent of the phenomenon of the oscillation's death due to the time- 
delay, although in this case the relevant dynamics is that of the averaged 
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global variables and not that of the individual neurons. 

We have derived the mean field approximation for the delayed coupled 
noisy system using the example of FitzHugh-Nagumo neurons in the excitable 
regime. It is expected that the approximations are equally valid for noisy 
delayed coupled type I excitable systems hke the Terman-Wang neurons, or 
for bursting neurons like the Hidmarch-Rose model. Also the approximation 
should be applicable under the same assumptions for neurons interacting by 
delayed chemical rather then electrical coupling. 

Acknowledgements This work is partly supported by the Serbian Min- 
istry of Science contract No. 141003. We should hke to acknowledge useful 
comments of the two referees. 

5 Appendix 

The system of equations (1) and (2) can be written in the form: 

edxi — {x — /3 — y + I)dt + c{< x{t — t) > —Xi)dt 
dyi = {x + b)dt + V2DdWi 

where: 

< x(t - r) >= - Y^i^iit - r) 

i 

The bracket < a; > is always used to denote the average over the units 
of the local variable Xi, which is, by the mean field assumption, for large N 
approximately equal to the average over the assumed Gauss distribution of 
the corresponding local variable Xi. 

Next we introduce deviations from the mean field: 

rixiit) =< x{t) > -Xi{t), UyXt) =< y{t) > -Viit). 

Deviations will always appear averaged over A^ i.e. in the form of < n^; > 
and < > so that the index i is in fact redundant. Correlations between 
centered moments are defined as 

Sx{t) =< nl >, Sy{t) =< nl >, u{t) =< n^Hy > 

Our goal is to derive the equations governing the evolution of the averages: 
X —< X >,Y —< y >, Sx, Sy,u. Due to the mean field assumption, these 
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averages can be computed as averages over the stochastic distributions of 
the local quantities, which are by assumption Gaussian. The equation for 
the derivatives of X =< x >,Y =< y >,Sx,Sy,u will contain averages of 
monomials in local variables of various orders. In order to handel these we 
shall need to use the formulas for the cumulant expansions up the fourth 
order the local quantities. The general formulas for the cumulant expansion 
can be found for example in [20] . Due to the assumed Gaussian distribution 
the third and the fourth order cumulants (and all of the higher order) are 
equal to zero, which will be used to express averages of monomial in local 
variables that appear in the evolution equations. 

Using the cumulant formulas one computes the following expressions 
which will be used to obtain the evolution equations: 

From the cumulant << x^y »= follows 

< x'^^yi >= Ys^ + YX^ + 2Xu. 
From the cumulant << x^y >>= follows 

< x^^yi >= 3s + 3X\ + YX^ + 3XYs^. 
Similarly one obtains: 

< X- > = + 

<x^i > = X^ + 3Xs^, 

<xt> = X^ + 6Xh, + 3sl, 

<Xiyi> = U + XY. (15) 

These expressions provide the necessary ingredients to obtain the equations 
(6). 

Takeing the average of the equations for x and y gives the first two equa- 
tions of the system (6). Next consider the equation for s^- 



= 2< X{t)X{t) - X{t)xi{t) - Xi{t)X{t) + Xi{t)xi{t) > 

^'^\x{t) - X{tf/3 - X{t)s,{t) - Y{t) + c{X{t - r) - X{t)] 



e 

+ - < Xi{ty - Xiitf/3 - Xi{t)y,{t) + cx,{t)X{t - r) - cx,{tf > 
e 

= '^[-X^{t)s,it) + s,it)-slit)-uit)-cs,{y)] (16) 
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which is the third equation (6). In the last equahty we used the expressions 
obtained from the cumulant formulas. 
Equation for Sy is obtained as follows: 

Sy^d< Y{tf - 2Y{t)yi{t) + yi{tf > /dt = -2Y{t)Y{t) -d< yi{t) > /dt. 

Using the Ito chain rule this becomes: 

- 2Y{t)[X{t) + b]+ < 2y(t)dyi{t)/dt + 2D> 

= -2Y(t)X{t) - 2Y(t)b+ < 2yi(t)xi(t) + 2yi{t)h + 2yi{t)V2DdWi + 2D> 

= -2Y{t)X{t) - 2Y{t)b + 2u{t) + 2X{t)Y{t) + 2Y{t)h + 2D 

= 2u{t) + 2D, (17) 

which is the forth equation (6). 

Similar calculations result in the u equation (6). 

u{t) = d<X{t)Y{t)-X{t)y,{t)-Y{t)x,{t)+x,{t)y,{t)> /dt 

= -X{t)X{t) - Y{t)X{t)+ < yi{t)xi{t) > + < yi{t)xi{t) >= . . . 

= ^u{t)[l-X^{t)-s^{t)-c]-^Sy{t) + s^{t). (18) 

References 

[1] Izhikevich E M 2005 Dynamical Systems in Neuroscience: The Geometry 
of Excitability and Bursting. The MIT Press. 

[2] Buric N and Todorovic D 2003 Dynamics of FitzHugh-Nagumo excitable 
systems with delayed coupling. Phys.Rev.E, 67 066222. 

[3] Mainen Z F and Scjnowski T J 1995 Reliability of spike timing in neo- 
cortical neurons. Science, 268 1503. 

[4] Haken H 2006 Brain Dynamics: Synchronization and Activity Patterns 
in Pulse-Coupled Neural Nets with Delays and Noise. Springer- Verlag, 
Berlin. 

[5] Linder B Garcia-Ojalvo J Neiman A and Schimansky-Geier L. 2004 Ef- 
fects of noise in excitable systems. Phys. Rep. 392 321. 



13 



[6] Buric N, Grozdanovic I and Vasovic N 2005 Type I vs Type II excitable 
systems with delayed coupling. Chaos, Solitons and Fractals, 23 1221. 

[7] Dhamala M, Jirsa V K Ding M 2004 Enhancement of Neural Synchrony 
by Time Delay. Phys. Rev. Lett. 92 074104. 

[8] Huang H and Cao J 2007; Exponential stability analysis of uncertain 
stochastic neural networks with multiple delays. Nonlinear Analysis: 
Real World Apphcations, 8 646. 

[9] Zhou Q Wan L and Sun J 2007 Exponential stabihty of reactiondiffusion 
generalized CohenGrossberg neural networks with time-varying delays. 
Chaos, Solitons and Fractals, 32 1713. 

[10] Buric N, Todorovic K and Vasovic N 2008 Synchronization of bursting 
neurons with delayed chemical synapses. Phys. Rev. E, 78 036211. 

[11] Janson N B Balanov A G and SchoU E 2004 Delayed feedback as a 
means of control of noise-induced motion. Phys.Rev.Lett. 93 010601. 

[12] Rodriguez R and Tuckwell H C 1996 Statistical properties of stochastic 
nonlinear dynamical models of single spiking neurons and neural net- 
works. Phys.Rev.E. 54 5585. 

[13] Tanabe S and Pakdaman K 2001 Dynamics of moments of FitzHugh- 
Nagumo neuronal models and stochastic bifurcations. Phys. Rev. E 63 
031911. 

[14] Zaks M A Sailer X Schimansky-Galer L and Neiman A B 2005 Noise 
induced complexity: from subthreshold oscillations to spiking in coupled 
excitable systems. Chaos. 15 026117. 

[15] Hasegawa H 2003 Dynamical mean-field theory of noisy spiking neuron 
ensembles: Apphcation to the Hodgkin-Huxley model. Phys. Rev. E. 68 
041909. 

[16] Hasegawa H 2004 Augmented moment method for stochastic ensembles 
with delayed couplings. II. FitzHugh-Nagumo model. Phys. Rev. E. 70 
021912. 

[17] Hasegawa H 2004 Augmented moment method for stochastic ensembles 
with delayed couplings. I. Langevin model. Phys. Rev. E. 70 021911. 



14 



[18] Arnold L 1998 Random Dynamical Systems, Berlin, Springer- Verlag. 

[19] Hale J and Lunel S V 1993 Introduction to Functional Differential Equa- 
tions, New- York, Springer- Verlag. 

[20] Gardiner C W 1985 Handbook of Stochastic Methods for Physics, Chem- 
istry and the Natural Sciences, Berlin, Springer- Verlag. 



15 



